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Abstract 

^ I This paper reports a method allowing for the determination of the 

shape of deformed fiber-like objects. Compared to existing methods, 
it provides analytical results including the local slope and curvature 
^ . which are of first importance, for instance, in beam mechanics. The 

' presented VIC (Virtual Image Correlation) method consists in looking 

CO i for the best correlation between the image of the fiber-like object and 

' a virtual beam image, using an algorithm close to the Digital Image 

Correlation method developed in experimental solid mechanics. The 
computation only involves the part of the image in the vicinity of the 
fiber: the method is thus insensitive to the picture background and the 
computational cost remains low. Two examples are reported: the first 
proves the precision of the method, the second its ability to identify a 
complex shape with multiple loops. 

^ • keywords : analytical shape, curve detection, virtual image correlation, 

^ fiber, filament, loop 



1 Introduction 

The determination of the shape and position of elongated objects such as 
hair [1], pulp fibers [2], needles [3], biological filaments [4, 5] or abiological 
objects [6, 7] is of interest in various fields of research. These objects may be 
bent by internal or external forces induced, for instance, by their own weight, 
by flowing fluids or by the interaction with solid surfaces. Image processing 
offers a set of techniques and algorithms that may be used for the detection 
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of such features. The basic technique consists in thresholding the image 
and skeletonizing it. An other possibility is to follow the line of maximal 
intensity (called "ridge") using the eigenvector corresponding to the highest 
eigenvalue of the Hessian of the (smoothed) image [8] or by using the minimal 
path method [9, 10]. Among the various families of techniques, a class of 
method is based on the Hough or the Radon Transform [11, 12, 13] which 
consist in transforming the images in such a way that segments become 
points, easily detected by image thresholding. Finally, the level set method 
also provides an efficient measurement of contours [14]. Yet, while these 
methods allow one to determine the shape of a curvilinear object, they do 
not allow one to estimate precisely its local curvature: this is however a 
key parameter which may be directly related to its mechanical state by the 
beam theory [15]. 

In the present method, the fiber shape is given by the optimal correla- 
tion between its physical image and a virtual beam. The latter has a mean 
line defined from a series expansion and a gray level which smoothly de- 
creases from the mean line to the borders. The correlation only involves 
the definition domain of the virtual beam (close to the physical one) that 
represents in general a much smaller number of pixels than the complete 
image. The correlation algorithm uses recent developments of the Digital 
Image Correlation techniques (DIG) and its application to mechanics [16]. 
This operation does not require any light intensity thresholding. 

The characteristics of the virtual beam are introduced in Sec. 2. The 
mathematical technique used for correlating the virtual beam onto the raw 
experimental image is proposed in Sec. 3. Sec. 4 details the determination of 
the initial parameters required by the method. Two practical applications 
of the Virtual Image Correlation (VIC) technique are reported in Sec. 5.1 
and 5.2. The first one consists of a straight cantilever elastic beam bending 
under its own weight: it is shown that the measurement is consistent with 
the theoretical result given by the beam theory. The second uses an experi- 
mental low resolution and noisy picture of a fiber transported (and curved) 
by a fiow in a fracture. The complex shape is recovered, demonstrating 
the robustness of the method with regards to loops, noise, luminance and 
contrast variations. Further developments of the VIC method are finally 
discussed in Sec. 6. 
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2 Parameterization of the virtual image 



The virtual image G consists of the virtual beam (in the sense of a plane 
curvilinear object in mechanics) with a length L and a width 2R. Any point 
X of the beam is parameterized by its curvilinear abscissa s G [0, L] along 
the mean line and its transverse distance r € [—R,R\ from it (see Fig. 1); 
the image G is not defined outside the beam's definition domain Dg- The 



local curvature is 7(5), the angle is 0{s), with 6*0 



0). The cartesian 




Figure 1: The virtual beam image and its coordinates system 



reference frame is (ei,e2), the tangent and normal vectors are respectively 
(t, v) and a point of the mean line is referred to as x (Eq. 1 and 2). There 
is no overlap, i.e. the points X are uniquely defined, since the local radius 
of curvature I/I7I is greater than the beam thickness R. 

X = x + ri/ (1) 

X = X0+ f T{i)di (2) 
Jo 

In order to give a finite dimension to the problem, the curvature 7(5) is 
described by a truncated series: 

TV 

7(s) = J^A„7„(5), (3) 

n=0 

where are the coefficients of the series, its order, 7„,(s) the dimension- 
less basis functions and s = s/ L ^ [0,1] the reduced curvilinear abscissa. 
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By integration, the angles are given by: 



N 



(4) 



n=0 



in which 




■s 



7n,(e)de. 



(5) 



Both 7„ and 9n depend only upon the choice of the series. They are esti- 
mated and stored prior to any computation. The set of parameters which 
must be determined are the magnitudes An and (xq,6o). 

The value of the gray levels g(X.) in the virtual image is chosen to be 
similar to that in the physical objects; it consists of a symmetrical function 
of the radius ^(X) = l{r) that continuously decreases from the mean line to 
the border: 



The present method does not require that the virtual image be closely similar 
to the physical one. Both the width and the gray levels of the virtual beam 
may be very different from the physical ones; this will be illustrated by the 
examples. 

3 Adjustment of the virtual beam to the physical 
picture 

As shown in the previous section, for a given set of function 7„, the shape 
of the virtual image G is fully described by the set of parameters V = 
{xq,9o, An} (a pseudo- vector of dimension + 4 whose components are 
referred as hereafter in this section). In this section, we describe the 
method developed to find the optimal set V-^, such that the beam image G 
falls best onto the image of the physical object that is contained in image F. 
Following [16], this is achieved by the minimization of the function <l>(Vfc): 



in which / is the luminance of the physical image F and dS = {l — jr)drds is 
the surface element. The domain Dg is fully contained inside F so that / is 
defined at any points. The optimal set V-^ corresponds to the minimum of ^ 
then the condition d$ = {d^/dVk)dVk = must be fulfilled for any variation 




(6) 




(7) 
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dVfc around Vjf . Using Eq. (7) together with dg/dVk = grad{g).dX./dVk, 
this condition leads to: 

where dDg is the boundary of the domain Dg i.e. the external boundary 
of the virtual beam, dl a differential line element and n a vector normal to 
the boundary and pointing outwards. Supposing that the iterative process 
is close to the solution, and that R is slightly greater than the width of 
the physical object, this boundary dDg is located in the background of the 
physical image (we neglect the boundary sides at s = and s = L as 
R « L). Then, assuming that the background is uniform, f\dDg (which 
represents the value of / along the boundary dDg) is constant. Furthermore, 
with the retained definition of g (Eq. 6), g\dDg = 0. Using the divergence 
theorem, we have: 

As div(X) = 2, the surface integral in Eq. (9) is constant (as 5 = 2RL) 
and its derivative with respect to equals zero. Thus, assuming hence- 
forth that the virtual beam boundary encloses the physical one and that the 
background of / is uniform, Eq. 8 reduces to: 

//_^(/-.)(grad„).||)d5 = 0. ,10) 

The next step consists in considering the Taylor expansion of g up to the 
first order: 

g{Vk + AVk) = g{Vk) + grad(g).— AUp (11) 
which, when introduced in Eq. (10), gives: 

This can be written as a matrix equation: 

MkpAVp = Lk, (13) 
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This equation represents a simple linear square matrix problem {N + 4 di- 
mension); its solution AVp is used to update the shape of the virtual beam. 
The iterative process is repeated until the value of $ decreases by less than 
a prescribed amount (10~^) between two steps. 

The term grad{g) .d^/ dV^ is involved in both the expressions of M^p 
and Lk- From Eq. (6) and the beam geometry follows: 



grad(5) = l'{r)v, 



(14) 



and, from Eq. (1): 



9X (9x 
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dVk dv, 'dYk-- ^^^^ 

The second term, collinear to r, does not need to be computed as it is 
orthogonal to grad(5'). The derivatives of dx/dV^ are obtained using Eq. (2) 
and (4): 



9x 

9X0,1 



ei, 



dx 



AC ^ 
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(16) 



The displacement of the virtual beam between two steps is close to ((9X/9Vfc)AVfc 



A 



V 
V 



Figure 2: Examples of unitary displacement fields. From left to right: 
dx/dxo^2 (vertical translation), Ox/OOq (rotation) and dx/dAo (uniform in- 
crease of the curvature) . 



where the fields d^/dVk represent unitary kinematic fields. Figure 2 shows 
an example of such fields (for clarity the fields are only represented at the 
mean line location, i.e. dx/dV^). They play the same role as the unitary 
displacement fields used in the DIG method [16, 17, 18]. 

The virtual image G is naturally discretized over the curvilinear frame 
(s,r) that does not correspond to the square grid of image F. To avoid any 
loss of information, the mesh size of the virtual image is much smaller. The 
computation of the second member of Eq. (12) requires one to project the 
luminance field / onto the mesh of G using, here, a cubic interpolation. 
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The overall length L of the beam is currently not straightforwardly ob- 
tained by the VIC method. The end of the fiber is detected from the sharp 
variation of d^/ds which occurs at the point where the virtual beam exceeds 
the end of the fiber. 



4 Determination of the initial conditions 

It has been shown that the VIC method, detailed in the previous section, 
requires that the boundary of the virtual beam initially surrounds the image 
of the physical fiber. This section briefly describes the method used to 
determine these initial conditions. 

The beam is first approximated by a sequence of straight segments. Their 
positions are given by a simplified VIC method in which the kinematic field 
is a single rotation around the first point of the segment. Finally the approx- 
imate shape is defined upon a collection of equally spaced points Xg (with 
^ q ^ Q). This gives the segments angles 9{sq) where Sq are the curvilin- 
ear abscissae of points Xg. These data are used to define the set of initial 
parameters = {-x.o,0o, An} that will be used in the first step of the VIC 
computation. The first term, xq, corresponds to the user-defined position of 
the initial point. The other terms are obtained by setting A-i = 6q/L and 
0-1 (s) = 1 in Eq. (4) which becomes: 

^ = f: ams). (17) 

n=— 1 

In particular, this equation applies to each value 9{sq) obtained from the 
previous approximative analysis, leading to: 

f^^EiA(l), as) 

n=—l 

where L = sg+i. From Eq. 5, the functions On depend only on the series 
functions 7^. Setting the order N of the series such asN + 2 = Q + l (the 
number of points x^), this equation represents a linear square matrix system 
whose resolution gives the terms An (where A-i = 9q/L). This defines 
and the VIC computation may start at this order (Q — 1), or at a lower one 
if is truncated. 
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5 Example of applications of the VIC method 



This section describes two examples especially chosen to illustrate, from a 
practical point of view, the different steps involved in the technique. The 
first example (Sec. 5.1) validates the accuracy of the method for a simple 
geometry; the second (Sec. 5.2) demonstrates its robustness in the case of a 
curled shape and a low quality image. 

5.1 The cantilever beam 

We proceeded to a simple experiment of a cantilever (clamped at one end) 
straight bar, bent under its own weight (Fig. 3a). The 2017-T4 aluminium 
bar has a length of 2459 mm and a radius of 4.95 mm. It is left free to bend 
under its own weight with one of its ends clamped in an horizontal chuck 
(Fig. 3b). A black curtain was placed in the back of the device in order 
to obtain a good intensity contrast and an uniform background (the latter 
condition is detailed in Sec. 3). A single light has been positioned close to 
the camera in order to illuminate equally both sides of the rod (in order 
to have a symmetrical luminance function /(r)). Once the bar is at rest 
(this may take a few minutes) an high resolution picture is captured using 
a Nikon D300 digital camera. This image (cropped from the original one) 
has 3897 x 1841 pixels. The field of view represents 2.33 x 1.09 m^ in the 
focal plane. 

For this image, we used Legendre series for the functions 7^ in Eq. (3). 
The Fourier series used in Sec. 5.2 provides similar results, but an higher 
order is needed. The Legendre series, in their shifted version {s G [0..1]), 
can be written as: 



in which the parenthesis refers to combination formulas. Yet, for our particu- 
lar example, because of the small contrast between the chuck maintaining the 
bar and the bar itself (see Fig. 3b), the abscissa xq^i was fixed so as to avoid 
any unwanted inclusion of the clamping device onto the virtual beam. The 
set of shape parameters (introduced in Sec. 3) reduces to = {xo,2) So,An}. 

Fig. 3 shows the result of the identification for an order = 3: the mean 
line of the virtual beam perfectly follows the middle of the aluminium bar 



7n = PnkS 



(19) 



where /c G [0, A^] and with [19]: 




(20) 
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1(d) 

Figure 3: (a) Aluminium bar bending under its own weight and identified 
mean line (white), (b-d) Magnified views of regions (b) close to the chuck, 
(c) midway along the bar and (d) at the tip. Solid black line: mean line 
obtained by the VIC method for a Legendre series of order 3. 
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along its full length. Due to its structure, the VIC method is not influenced 
by large and illuminated objects present in the foreground of the picture (like 
the stones on the bottom left of Fig. 3a). The virtual beam diameter was 
set at 20 pixels (the physical bar diameter is around 10 pixels). The virtual 
beam mesh has 61 x 8634 points (approximatively three times thinner than 
the image grid). The total time required for the computation was about 3 
seconds per iteration using a 2.8 GHz dual-core computer and an algorithm 
implemented on Matlab. 

We now compare these results to the prediction of the elastic beam theory 
[15]. According to it, the flexural moment M{s), given by static equilibrium, 
is proportional to the local curvature: 



where the abscissa xi is measured along the horizontal axis. This system 
is solved using a numerical iterative method and the physical properties 
p = 2700 g/cm^ and E = 72 GPa [20]. The Young modulus was confirmed 
by a three-point bending test performed on the actual specimen, leading to 
72.6 GPa. 

Fig. 4a shows that the results given by the beam theory match very well 
the measurements from the VIC method, in term of displacement in the 
xi (horizontal) X2 (vertical) plane. Fig. 4b presents a magnified view of 
the same result and only small differences are observed: the mean quadratic 
vertical distance between the two solutions is approximatively of 2 pixels (i. e. 
1.2 mm). A similar discrepancy has been observed in a second experiment in 
which the chuck was rotated by an half turn. Figs. 4c and 4d demonstrate 
that the VIC method determines correctly the slopes and curvatures. In 
particular, both the initial slope 9{s = 0) = and the final curvature 7(s = 
L) = are correctly determined (note that their value were not imposed by 
the Legendre polynomial sequence used). 

The choice of the order N results from a classical compromise between 
simplicity and accuracy. An order too low (here = 1) may lead to a loss of 
correlation between F and G (the boundary dDg does not enclose anymore 
the image of the object). Increasing the order monotonically decreases the 
value of the correlation function (Fig. 5). When the order becomes too high 
(for A^ > 12 in this example) the matrix M^p in Eq. (13) becomes ill-defined 
while the level of <I> asymptotically reaches a minimum value. 



M{s) 
6(0) 




(21) 



(22) 
(23) 
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Figure 4: Comparison between the beam theory (sohd hne) and the VIC 
method (circles), (a) Displacement, (b) Discrepancy between the two deter- 
minations of the vertical displacement, (c) Angle (radians), (d) Curvatures. 
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Figure 5: Dependence of the final value of the correlation function ^ on the 
order N. 
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Figure 6: a) and b) : the physical beam image in the unwrapped coordinates 
of the virtual beam for = 3 and = 8. 



A visual comparison of the results obtained for two different values of 
A^ is displayed in Figs. 6a and 6b. They represent the physical image in 
the reference frame (s, r) of the virtual beam (represented unwrapped for 
clarity). In the case of an ideal correlation, whatever the shape of the fiber, 
this figure should be symmetrical with respect to the (straight) mean line. 
Because of the magnification along the r axis, one can see that the image in 
Fig. 6a obtained for a low order A^ = 3 is slightly wavy while the image in 
Fig. 6b obtained for A^ = 8 is perfectly symmetrical. The inhomogeneities 
of the illumination and of the background are visible but do not influence 
the final result. 

5.2 The fiber transported by a fiuid flow in a fracture 

This second example illustrates the robustness of the technique with regard 
to the quality of the physical image F and to the complex shape of the 
object. 

This example deals with the transport of a fiber by a fiow fiuid in a 
transparent fracture made of two rough surfaces [21] (Fig. 7). Because of 
the fracture roughness, the free space has a complex geometry leading to 
a disordered flow velocity field. As a result, when a fiber of diameter (0.3 
mm) close to the fracture aperture (0.75 mm) is inserted into the fracture, 
it is subject to spatially variable forces from the fluid resulting in strong 
deformations of the thread. Fiber transport is only possible at large fluid 
velocities so that images need to be captured at high rates with a poor light 
intensity contrast. 

Several attempts using thresholding operations and spline interpolations 
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20 40 60 80 1 00 1 20 1 40 1 60 1 80 200 220 (d) 



Figure 7: Fiber transported by a fluid flow in a fracture (units are in pixels), 
(a) Physical image, (b) Same image and identified mean line (white line), 
(c) Detailed view: identified mean line (white line) and border of the virtual 
image (white dots), (d) Virtual beam. 
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were made to extract the center line of the fiber but none of them gave a 
solution well fitted to the whole profile of the fiber at a given time. This 
was especially obvious in situations of large deformation with buckling. The 
noise level of the picture was a major problem for thresholding methods 
(the light intensity of spots in the background is close to that of the pixels 
corresponding to the fiber). 

The experimental pictures were them reanalyzed using the VIC method. 
A Fourier series containing = 100 terms was necessary to describe this 
complex shape (the previously used Legendre series lead to numerical prob- 
lems due to large values of Pnk in Eq- 20). Fig. 7c shows that the fiber 
diameter is close to one pixel then a slightly larger width, R = 1 pixel, has 
been selected for the virtual beam. Similarly to previous example, the ab- 
scissa xo,i of the first point (on the left of Fig. 7a) has been fixed in order 
to avoid unwanted edge effects. 

The Fig. 7 clearly shows that the VIC method gives an accurate estima- 
tion for the central line of the fiber. Despite the small radius of the object 
and the high noise level, the central line collapses everywhere onto the fiber. 
Wide meanders as well as small deviations from a straight line are perfectly 
followed. Moreover, even the loops of the fiber can be extracted by this 
method. 

6 Conclusion and discussion 

The Virtual Image Correlation method allows one to identify precisely the 
shape of a fiber from its image. The virtual beam fully encloses the physical 
image of the fiber so that all the information contained in it is used (and 
not only its brightest pixels) ; this feature contributes to the precision of the 
method. Furthermore, pixels outside of the definition domain are not used: 
this speeds up the computation and makes it insensitive to possible artefacts 
in the background. The method can deal with noisy images and strongly 
curved fibers; it can already be used in various domains, from biology to me- 
chanical engineering. The analytical identification of the curvatures given 
by the VIC provides informations on the mechanical state of the fiber: it 
is possible to identify its fiexural properties from the knowledge of external 
forces. Conversely, if the mechanical parameters are known, one can com- 
pute (using an inverse problem approach), the external forces acting on the 
fiber. 

From a theoretical point of view, it would be interesting to include the 
length of the fiber in the optimization process. The examples provided in 
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the present paper use either Fourier or Legendre series decompositions for 
the representation of the mean hne. Other series such as the Chebyshev 
polynomials or other functions used for Computed Aided Design (Bezier 
curves, splines...) may also be used. Following [22, 16], if the method is 
used for identifying the mechanical properties of the fiber, one should define 
the main term of the series from the theoretical beam equation. 

Time sequences will also be considered: the variation of images with time 
allow one to increase the precision of the analysis and to study the dynamic 
of the buckling process. Three-dimensional analysis is also envisaged and 
will use at least two pictures taken from different angles of view. 

The VIC method can also be applied to images of other curvilinear 
shapes (non fiber objects) in various scientific domains in which the cur- 
vature has a major role, for example thermal or diffusion fronts in fluid me- 
chanics, chemistry... It may also be useful in the "pure" image processing 
field in some case of line or pattern recognition. With a slight modification 
of the virtual beam luminance definition (a step shape), the method may 
also be applied to edge detection problems. 
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